home *** CD-ROM | disk | FTP | other *** search
/ El Mac 9 / El Mac 9.iso / Shareware / Applications / MathPad 2.4 / Examples / Navigation < prev    next >
Encoding:
Text File  |  1995-04-30  |  2.2 KB  |  81 lines  |  [TEXT/MPad]
  1. -- Great circle route from 'here' to 'there'
  2. -- by Hank.Dolben@UNH.edu 1994 May 18
  3. ---------------------------------------------
  4. include ":incl:rotations"
  5.  
  6. places = {Anchorage, BuenosAires, CapeTown, Honolulu, London,
  7.        Moscow, NewYork, Singapore, Sydney, Tokyo, Wellington}
  8.  
  9. here=London; there=Anchorage
  10.  
  11. az    := heading(here,there):;            az:-15.6
  12. alpha := angle(here,there):; distance(alpha):4457.0
  13.  
  14. -- For a nice background, use the world map picture from the system 7 distribution Scrapbook File. Paste it into the PLOT and comment out the following line.
  15. plot place(places)
  16.  
  17. -- locations are given in {+N latitude, +E longitude} degrees.
  18. Anchorage  ={ 61,-150}
  19. BuenosAires={-35,- 58}
  20. CapeTown   ={-34,  18}
  21. Honolulu   ={ 21,-157}
  22. London     ={ 52,   0}
  23. Moscow     ={ 56,  38}
  24. NewYork    ={ 41,- 74}
  25. Singapore  ={  1, 103}
  26. Sydney     ={-34, 151}
  27. Tokyo      ={ 36, 140}
  28. Wellington ={-41, 174}
  29.  
  30. plot place({here})
  31.  
  32. s := rotation(here[1],here[2],90):
  33. inc := multiply(rotation(incangle,0,az),rotation(0,0,-az)):
  34.  
  35. plot {step,c[1]}; plot place({there})
  36.  
  37. step = r:={1,0,0} when X=0,
  38.        r:=transform(inc,r), b:=transform(s,r), c:=geog(b),
  39.        c[2]
  40.  
  41. place(L)[i] = {L[i,2],L[i,1]}  -- swap for plot. X=lon Y=lat
  42.  
  43. increment=300; incangle=increment/LenPerDeg
  44. Xmin=-180;  Xmax=180;  Xsteps=trunc(alpha/incangle); Xdiv=45;
  45. Ymin= -89;  Ymax= 89; Ydiv=45;
  46.  
  47. -- the great circle heading from A to B
  48. -- is the azimuth of B-A at A
  49. heading(A,B) = azimuth(head(A,B))
  50.  
  51. -- the arc length of an angle A
  52. distance(A) = A*LenPerDeg
  53.  
  54. LenPerDeg = radius*π/180
  55.  
  56. -- the radius of earth (pick distance units)
  57. radius = 3960
  58.  
  59. -- azimuth is measured east of north
  60. -- in the spherical system V[2] is south, V[3] is east
  61. azimuth(V) = atan2(V[3],-V[2])
  62.  
  63. -- B-A at A
  64. head(A,B) = rotate(A[2],-A[1],-90,cart(B)-cart(A))
  65.  
  66. -- angle between locations A and B
  67. angle(A,B) = arccos(dot(cart(A),cart(B)))
  68.  
  69. -- cartesian coordinates for unit vector at A={lat,lon}
  70. cart(A) = {sin(90-A[1])*cos(A[2]),sin(90-A[1])*sin(A[2]),cos(90-A[1])}
  71.  
  72. -- {lat,lon} of cartesian unit vector A
  73. geog(A) = {90-arccos(A[3]),atan2(A[2],A[1])}
  74.  
  75. -- limit the argument of acos (fixes numerical error)
  76. arccos(cosine) = acos(limit1(cosine))
  77.  
  78. limit1(a) = a when abs(a) ≤ 1, sign(a) otherwise
  79.  
  80. sign(a) = -1 when a < 0, 1 otherwise
  81.